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ABSTRACT 


> 


The  problem  of  modeling  and  processing  radar  images  for  tracking  ap- 
plications has  been  considered.  Digital  simulation  experiments  have  been 
carried  out  for  area  correlation  of  a reference  target  image  with  on-board 
radar  acquired  images.  Effects  of  misregisterations  of  radar  height  (scaling) 
and  orientation  (rotation)  and  those  of  high  pass  filtering  of  images  have 
been  studied.  Use  of  image  models  based  on  finite  difference  approximation 
of  partial  differential  equations  (PDEs)  has  also  been  studied. 


KEY  WORDS:  Image  Modeling;  Area  Correlation;  Radar  Image  Processing 
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I.  SUMMARY 


During  the  duration  of  this  project  the  following  progress  has  been 

made. 

1.  The  problems  associated  with  the  overall  goal  of  Forward  Looking  Radar 
Image  Modeling  and  their  real  time  simulation  have  been  identified. 

The  basic  problems  are: 

a)  Modeling  and  Identification  of  the  Point  Spread  Function  (PSF) 
of  the  FLR  seeker. 

b)  Stochastic  Modeling  of  Image  Fields  and  their  use  in  image  re- 
storation and  system  identification 

c)  Hardware  encoding  of  FLR  image  model  for  real  time  RF  simulation 

2.  A mathematical  formulation  of  a typical  FLR  image  has  been  developed. 
The  FLR  image  can  be  synthesized  digitally  using  this  model  and  a 
reference  Side  Looking  Radar  (SLR)  image.  The  model  can  also  be 
identified  if  the  statistical  properties  of  SLR  image  are  known. 

The  underlying  identification  scheme  would  require  a statistical 
model  for  two  dimensional  image  fields.  The  FLR  images  are  used 
for  re-entry  vehicle  guidance  via  area  correlation  techniques.  For 
a high  signal  to  noise  ratio  of  correlation  peaks  the  FLR  image 
(which  is  generally  quite  blurred)  should  be  restored  (or  enhanced). 
Image  restoration  algorithms  have  been  studied  in  this  period  using 
some  new  image  models  and  some  new  algorithms  have  been  shown  to  be 
associated  with  these  models.  (see  attached  list  of  publications). 
Further  work  is  required  to  consider  restoration  schemes  which  in- 
clude space  varying  nature  of  the  point  spread  function  (PSF)  of  the 
FLR  imaging  system. 
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3. 


The  problem  of  stochastic  Image  modeling  has  been  studied  in  the 


framework  of  stochastic  partial  differential  equations.  Our  re- 
sults show  that  these  models  are  quite  effective  and  are  often 
better  than  conventional  image  covariance  models.  These  models 
have  been  applied  for  restoration  of  noisy  images.  Further  work 
is  required  for  application  of  these  models  to  restoration  of 
FLR  images  and  for  further  refinement  models. 

4.  Collaborative  efforts  were  made  with  the  Radio  Frequency  Simula- 
tion System  (RFSS)  laboratory  of  the  Advanced  Simulation  Center 
(ASC)  of  the  U.S.  Missile  Research  and  Development  Command 
(MlRADCOl)  for  hardware  in  the  loop  simulation  of  the  FLR  image 
data.  These  efforts  have  produced  a design  of  methodology  of  RFSS 
configuration  to  generate  RF  signals  which  represent  FLR  image 
data.  Harware  design  for  transmission  and  reception  of  these 
signals  still  needs  to  be  accomplished. 
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II.  FORWARD  LOOKING  RADAR  IMAGE  MODELING  AND  SIMULATION 


2. 1 Introduction 

It  is  well  known  that  for  missile  and  aircraft  guidance  the  area  cor- 
relation between  a target  reference  image,  already  stored  in  the  vehicle, 
and  a target  image  acquired  by  an  on-board  sensor  could  be  utilized.  The 
cross-correlation  function  of  the  two  images  is  first  obtained  and  the  peak 
value  of  chis  function,  indicating  a match  between  the  two  images  is 
located.  The  peak  location  is  used  to  guide  the  vehicle  in  the  direction 
of  the  target.  It  is  assumed  that  a significant  area  of  the  target  image 
is  within  the  range  of  the  on-board  reference  image;  otherwise,  no  peak 
(or  a false  peak)  could  be  located.  Typically,  this  technique  is  used  in 
the  terminal  phase  of  the  flight. 

For  many  terminal  guidance  applications  the  reference  image  is  often 
generated  by  a side  looking  radar  (SLR)  or  a high  resolution  sensor,  while 
the  on-board  sensor  is  a forward  looking  radar  (FLR)  with  significantly  re- 
duced resolution.  Application  of  area  correlation  for  guidance,  therefore, 
requires  an  understanding  of  interactions  between  the  actual  scene  and  the 
radar  sensing  process.  Also  important  for  obtaining  a good  correlation 
peak  is  the  proper  regls teration  of  the  two  images.  In  practice  two  forms 
of  misregisterations , namely,  scaling  and  rotation  of  one  image  with  respect 
to  another  are  common.  So  it  becomes  Important  to  know  the  effects  of 
these  misregisterations  on  the  detection  and  the  location  of  the  peak.  Here 
we  report  the  results  of  the  computer  simulation  experiments  carried  out  to 
study  the  effects  of  FLR  parameters,  seal ing, rotation  and  the  filtering  of 
FLR  images  on  the  area  correlation. 

In  section  2.2  we  describe  the  theoretical  formulation  of  the  experi- 
ments and  the  necessary  assumptions.  Modeling  of  FLR  images  from  a given 
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SLR  image  is  briefly  described  in  Section  2.3.  Experimental  results  of 
computer  simulations  are  summarized  in  Section  2.4.  Conclusions  and  sug- 
gested future  work  are  given  in  Section  2.5. 


2. 2 Area  Correlation 

Let  u(x,y)  and  v(x,y)  denote  the  reference  and  the  on-board  generated 
images  respectively.  Then  their  area  correlation  function,  denoted  by  r(x,y), 
is  given  by 

00  00 

r(x,y)  = ''  V u(x',y')v(x  + x , y + y')dx'dy‘'.  (1) 

'-’-CO‘^-00 

In  practice  the  images  are  available  only  over  a finite  area  and  hence  the 
integration  in  (1)  is  performed  only  on  the  available  area.  For  digital 
processing  we  approximate  the  integration  in  (1)  by  a summation  given  by 

” ^ N N 

r(m,n)  =2  11  u(k,ji)v(k+m-l,  .^+n-l),  -r  - 1 5 m,n  ^ — , (2) 

k=U=l  ^ 

where  the  available  reference  SLR  signal  has  been  digitized  to  N x N samples. 
Evaluation  of  r(m,n)  using  (2)  requires  the  knowledge  of  one  of  the  images 
over  a larger  area  than  the  other.  If  both  images  are  available  over  the 
same  area,  as  is  assumed  in  our  experiments,  then  replacing  the  function 
v(i,j)  by  its  two  dimensional  periodic  repetition,  allows  one  to  compute 
(2)  efficiently  using  discrete  Fourier  transform.  If  U(m,n)  and  V(m,n) 
denote  the  discrete  Fourier  transform  (DFT)  of  sequences  u(k,£)  and  v(k,X) 
respectively,  then 

R(m,n)  = U*(m, n)V(m,n)  , I S m,n  ^ N (3) 

gives  the  DFT  of  r(k,ji),  where  * dentoes  the  complex  conjugate.  The  dis- 
crete correlation  function  r(k,£)  could  then  simply  be  computed  by  taking 
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the  Inverse  DFT  of  R(m,n).  The  use  of  fast  Fourier  transform  (FFT)  algorithm 

4 2 

reduces  the  computation  of  (2)  from  N to  0(N  log2N).  It  should  be  noted 
that  r(k,f)  computed  as  above  is  not  centered  around  the  origin.  This 
could  easily  be  done  by  shifting  the  samples  in  each  coordinate  axes  by 
^ and  re  indexing. 

2.2.1  Effect  of  Scaling 

The  effect  of  misregisteration  of  radar  heights,  at  which  the  reference 
and  the  on-board  images  are  generated,  is  equivalent  to  scaling  the  axes  of 
one  of  the  images  relative  to  the  others.  We  assume  the  scaling  of  the 
reference  image.  A digital  scaled  image,  with  a scale  factor  of  s > 0, 
could  be  approximated  by 

Ug(k,£)  = u(LskJ  , Ls^J)  (4) 

where  [aj  denotes  the  highest  integer  less  than  or  equal  to  a. 

2.2.2  Effect  of  Rotation 

Another  common  form  of  misregisteration  is  the  rotation  of  the  live 
image  with  respect  to  the  reference  image  along  the  vertical  axis  of  FLR. 

As  FLR  images  are  radially  scanned  the  rotation  could  be  easily  accomplished 
in  polar  coordinates.  Let 

v(r,9)  = v(rcos6,  rsin0) 

then  an  FLR  image  rotated  by  an  angle  cp  with  respect  to  the  reference  image 
in  clockwise  direction  is  given  by 

V (r,e)  = v(r,e-cp) 


and 
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/ 2 2 

v^(x,y)  = v^(V^x  +y  , arctan(y/x)). 

2.2.3  Enhancement  of  FLR  Images 

As  noted  earlier,  the  on-board  generated  FLR  Images  have  a poor  reso- 
lution compared  with  the  reference  images.  It  is  intuitively  clear  then, 
and  also  would  be  seen  in  the  next  section,  that  the  resulting  correlation 
function  would  not  yield  a desirable  sharp  peak  and  could  even  give  a false 
peak.  To  remedy  this  it  is  desirable  to  enhance  by  some  sort  of  high  pass 
filtering  of  the  FLR  images  before  correlating  them  with  reference  images. 
There  are  various  image  enhancement  techniques  that  are  available  in  the 
literature.  A class  of  techniques,  known  to  be  effective  for  blurred  images 
performs  filtering  in  the  Fourier  domain.  This  could  be  profitably  exploited 
in  our  case  as  we  are  already  taking  the  Fourier  transform  to  get  the  cor- 
relation function.  If  V(m,n)  is  the  DFT  of  the  discrete  FLR  image,  then  a 
typical  Fourier  domain  filter,  is 

Vj.(m,n)  = F(V(m,n)) 

where  F denotes  a filter  function,  which  could  be  linear  or  non-linear.  An 
example  of  linear  filter  is  the  Wiener  filter.  A nonlinear  high  pass  filter 
is,  for  example,  the  so-called  cx-f liter  given  by 

V^(m,n)  = |v(m,n)!"  eJ2narg(V(m,n))  q < ^ 1 (5) 

and  (3)  becomes 

R^(m,n)  = U*(m,n)V^(m,n)  . 

2.2.4  Performance  Evaluation 

In  practice  v(x,y)  and  u(x,y)  are  not  zero  mean  images.  If  the  mean  is 
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comparable  to  the  variance  of  the  Image  It  adds  a significant  constant  bias 
to  the  correlation  function  and  the  relative  sharpness  of  the  peak  is  de- 
creased. So  it  is  desirable  to  make  one  of  the  Images  zero  mean,  which 
in  our  case  could  be  easily  accomplished  by  setting  the  first  element  of  the 
DFi’  of  one  of  the  images  to  zero.  To  study  the  effects  of  FLR  parameters, 
misregisterat ions, and  a-filter  on  the  area  correlation  we  define  a 
criterion,  called  signal  to  noise  ratio  (or  S/N),  as 

„ , _ Value  of  the  Positive  Peak 

Standard  Deviation  of  the  Correlation  function 

2. 3 Modeling  of  FLR  Images 

For  analyzing  area  correlation  techniques  it  is  necessary  to  have  a 
model  of  the  PSF  of  the  FLR  image.  The  model  used  in  our  experiments  is 
based  on  the  theory  (which  is  also  supported  by  practical  data)  that  terrain 
reflectivity  is  independent  of  the  viewing  aspect  angle  except  for  extreme 
angles. 

Figure  1 shows  the  FLR  geometry.  Scanning  an  offset  antenna  pattern 
around  the  vertical  axis  results  in  a ring  or  doughnut  shaped  coverage  of 
ground.  At  any  scan  position,  the  resolution  cells  are  bounded  by  the  an-  j 

tenna  half-power  elliptical  contour  and  a radial  width  equal  to  one-half 
of  the  vertical  ground  intercept  of  the  radar  pulse  length.  This  width  is 
a secant  function  of  the  range  vector  depression  angle  0.  The  parameters 
of  the  half-power  ellipse  are  also  given  in  Figure  1.  In  polar  coordinates 
the  FLR  signal  could  be  expressed  in  terms  of  the  reference  scene,  u(r,cp)  by 

v(r,cp)  = \ \ u(r+s,  cp  •+^)dYds  (6) 

^-cp^(r)  '’-i^(r) 
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Antenna  half 
power  ellipse 


Cross  Section  of  Radar  Beam 


where  r = htanO,  l-s  the  angular  width  of  the  elliptical  contour  from 

Its  major  axis,  ij^(r)  and  correspond  to  the  Inner  and  outer  ground 

Intercept  of  radar  pulse  width  around  the  point  and  h Is  the  height  of  the 


radar. 


Equation  (6)  could  be  generalized  by  associating  a non-unity  kernel 
(or  weighting  function)  as 


v(r,cp)  = \ \ w(r,s;cp,Y)u(r+s,c(>+^)dYds 

t_>  /^\tJ  § /^\ 


-cp^(r)-’-A^(r) 


Euqatlon  (7)  could  be  simplified  by  making  the  assumptions  that  the 
weighting  function  is  separable  and  is  Invariant  with  respect  to  cp,  giving 


(r,cp)  = \ \ w (r,s)w  (Y)u(r+s,cp45i:)dYd£ 


-cp^(r) 

Note  that  the  point  spread  function  (PSF)  as  defined  above  Is  spatially  vary- 
ing so  that  conventional  Fourier  domain  Wiener  filtering  is  not  applicable. 

From  the  statistics  of  reference  image  and  FLR  Image  the  point  spread 
function  could  be  identified  by  assuming  some  models  of  the  weighting  func- 
tions, e.g.  Gaussain  density,  and  then  identifying  the  parameters  of  these 
models  by  using  mean  square  or  other  criteria- 


2.4  Experimental  Simulation  Results 


Figure  2 shows  a 256  x 256  remotely  sensed  six-bit  image.  Considering 
this  as  only  a single  quadrant  of  a reference  image, a 90*  segment  of  an 
FLR  image  was  simulated  using  Eqn.  (6).  The  blur-factor  as  reported  in 
our  experiment  is  roughly  a measure  of  average  number  of  samples  of  the 
reference  image  summed  to  get  the  intensity  at  each  sample  of  FLR.  Also, 
the  FLR  images  have  been  normalized  to  have  the  same  mean  square  energy  as 


figure  2 


Square  ERTS  Image  Used  as  SLR  Image 


FIGURE  3 


SLR  Image  of  Figure  2 With  Axes  Scaled  by  a Factor  0.98 


the  reference  image  within  the  area  covered  by  the  FLR  image. 
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Figure  3 shows  a scaled  reference  image  corresponding  to  Figure  2 with 
a scale  factor  of  0.98.  To  simulate  a doughnut -shaped  quadrant  of  an  FLR 
image  a doughnut-shaped  quadrant  of  the  square  reference  image  was  created 
and  sampled  in  polar  coordinates  with  4 samples  per  degree  in  azimuth  and  200 
samples  along  the  elevation.  Figure  4 shows  this  corresponding  to  Figure  2. 
From  this  three  different  FLR  images  as  discussed  in  Section  2.3  with  para- 
meters shown  in  Table  I were  obtained.  Figure  5 shows  one  such  image. 

Both  in  Figure  4 and  Figure  5 the  image  intensity  outside  the  doughnut  is 
replaced  by  its  mean  intensity.  For  calculating  the  correlation  function  the 
polar  FLR  images  were  transformed  to  cartesian  coordinates  square  images. 
Figure  6 shows  the  FLR  image  of  Figure  5 with  a 2 ® misregistration.  Figure 
7 and  Figure  8 show  the  logarithm  of  the  magnitude  of  the  discrete  Fourier 
transform  of  the  images  of  Figure  2 and  Figure  5 respectively.  Figure  9 
and  Figure  10  are  the  high  pass  filtered  images  corresponding  to  the  FLR 
image  of  Figure  5. 

The  results  of  digital  correlation  experiments  are  sumnarized  in  Table 
II.  Figure  11  and  Figure  12  show  the  correlation  function  and  its  cross- 
section  for  images  of  Figure  2 and  Figure  4.  To  improve  correlation  peaks 
the  reference  image  was  high  pass  filtered  (a  = .05). 

Figures  13-16  show  the  three  dimensional  plots  of  the  correlation 
functions  for  some  of  the  entries  of  Table  II.  Here  every  fourth  sample  of 
the  256  X 256  correlation  function,  has  been  plotted.  Figures  17-23  show 
the  horizontal  and  the  vertical  cross-sections  of  some  of  the  correlation 
functions  across  the  peak. 

Figure  24  shows  the  effect  of  FLR  beam  pulse  width  product  (or  equiva- 
lently blur  factor)  on  the  signal  to  noise  ratio  for  the  original  and  filtered 
reference  image.  Figure  25  shows  the  effect  of  scaling  the  axes  of  the 
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reference  image  on  the  signal  to  noise  ratio.  Figure  26  shows  the  effects 
of  rotating  the  FLR  image  and  the  combined  effect  of  scaling  and  rotation 
' on  the  signal  to  noise  ratio.  Finally,  Figure  27  shows  the  effect  of  the 

filter  parameter  a on  the  signal  to  noise  ratio. 

These  experiments  show  that  high-pass  filtering  improves  the  correlation 
results  considerably.  Fro".  Figure  24  we  see  that  the  signal  to  noise  ratio 
of  the  cross-correlation  between  the  reference  image  and  FLR  images  is  much 
lower  than  that  of  the  autocorrelation  of  reference  image  even  for  small 
amount  of  blur.  This  is  mostly  because  of  the  space-variant  nature  of  the 
FLR  blurring  mechanism  rather  than  the  low  resolution  of  the  FLR  images. 

Figure  25  shows  that  the  signal  to  noise  ratio  drops  sharply  with 
the  scaling  of  axes  and  is  almost  linear.  In  the  presence  of  large  rota- 
tional misregisteration  the  scaling  has  only  marginal  effect,  as  could  be 
seen  from  Figure  26. 

2. 5 Conclusions  and  Future  Work 

As  demonstrated  by  results  in  the  previous  section  the  enhancement  of 
the  on-board  generated  FLR  images  using  Fourier  domain  high  pass  filtering 
is  quite  useful  in  a more  unambiguous  location  of  the  correlation  peak. 
Although  this  simple  and  adhoc  technique  gives  good  results,  it  does  not 
perform  the  optimum  filtering  of  the  point  spread  function  of  the  FLR,  which 
is  space-variant  along  elevation.  A better  technique  is  to  identify  the 

models  of  the  PSF  and  then  do  easily  Implementable  optimum  or  sub-optimum 

I 

filtering. 

I 

j Stochastic  modeling  of  the  reference  image  and  the  FLR  image  is  another 

1 area  for  future  research.  A considerable  progress  has  been  made  recently  in 

I this  area  for  certain  class  6f  images  (see  Appendixes  I and  II).  These 
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models  could  be  used  for  Identification  of  PSF  and  also  In  applying  the 
statistical  detection  theory  to  identify  and  locate  the  correlation  peak. 
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FIGURE  4 


Ring-shaped  (Radial)  SLR  Image  obtained  from  Image  of  Fig.  2 
with  intensity  in  the  area  outside  the  ring  set  equal  to  the 
mean  intensity  of  the  image  inside  the  ring 
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FIGURE  7 

Logarithm  of  the  magnitude  of  the  Fourier  transform  of 
SLR  Image  of  Figure  2 


FIGURE  8 


Logarithm  of  the  magnitude  of  the  Fourier  transform  of 
FLR  Image  of  Figure  5 
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FIGURE  9 

cv-Filtered  image  of  the  FLR  Image  of  Figure  5,  a = 0.8 


FIGURE  10 

ff-Filtered  image  of  the  FLR  image  of  Figure  5,  o'  = 0.6 
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Horizontal  and  vertical  cross-aectlonal  plots  of  the  correlation 
function  shown  In  Fig.  H 


Figure  13:  Three  dimensional  plot  of  the  cross  correlation  function  of 
the  SLR  and  FLR-2  images 


Figure  15:  Three  dimensional  plot  of  the  cross  correlation  function  of 
the  scaled  (2%)  and  filtered  (o'  = .5)  SLR  image  and  rotated 
(2®)  FLR-1  image 
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Figure  17:  Horizontal  and  vertical  cross-sect  ions  of  the  3-D  plot 
of  Fig.  13  across  the  peak 
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Figure  19:  Horizontal  and  vertical  cross-sections  of  the  cross 

correlation  functlbui  of  rotated  (2®)  FLR-l  and  filtered 
(O'  = .5)  SLR  Images  across  the  peak 
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Figure  20;  Horizontal  and  vertical  cross-sections  of  the  cross 

correlation  function  of  FLR-l  and  scaled  (27J  and  filtered 
(o'  = .5)  SLR  Images  across  the  peak 


Figure  21:  Horizontal  and  vertical  cross-sections  of  the  cross 

correlation  function  of  FLR-2  and  scaled  (2%)  and  filtered 
(or  “ .5)  SLR  Images- across  the  peak 
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Figure  22:  Horizontal  and  vertical  cross-sections  fo  the  3-D  plot 
of  Figure  15  across  the  peak 
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Figure  23:  Horizontal  and  vertical  cross-sections  of  the  3-D 
Plot  of  Figure  16  across  the  peak 
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Figure  24:  Signal  to  noise  ratio  of  the  correlation  peak  vs.  the 
Blur  factor 
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Effect  of  scaling  on  the  signal  to  noise  ratio  of  the 
correlation  peak 
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III,  STOCHASTIC  IMAGE  MODELING  AND  APPLICATIONS 

The  simplest  stochastic  model  of  an  image,  at  least  conceptually,  is 
the  joint  probability  density  function.  However,  for  practical  reasons 
(of  dimensionality),  it  is  often  convenient  to  work  with  a two  dimensional 
stochastic  difference  (or  differential)  equation  representation.  State 
variable  models,  which  convert  two  dimensional  images  into  a sequence  of 
one  dimensional  vectors,  have  been  developed  and  implemented  for  a variety 
of  image  restbration  problems  by  Nahi,  Silverman,  et.  al  [1,2].  Ekstrom 
and  Woods  [3J  have  suggested  a two  dimensional  spectral  factorization  techni- 
que for  developing  non-symmetric  half  plane  (NSHP)  models  which  can  be  used 
for  recursive  image  processing  applications. 

We  have  demonstrated  that  the  theory  of  Partial  Differential  Equations 
(PDEs)  provides  a strong  framework  for  studying  many  problems  in  image  pro- 
cessing [4,5].  Many  of  these  models  have  been  found  to  be  superior  to  a 
widely  used  separable  covariance  model  for  images  in  image  restoration  ap- 
plications [61.  These  models  are  divided  into  three  classes  according  to 
the  classification  of  PDEs  as  hyperbolic,  parabolic  and  elliptic  equations. 

It  has  been  shown  in  [4]  that  these  classes  of  PDEs  could  realize  different 
types  of  correlation  functions.  Thus,  a given  correlation  function  may  be 
best  realized  (in  terms  of  the  minimal  order  of  the  PDE)  by  a particular 
class  of  PDEs,  For  example,  certain  isotropic  correlation  functions  which 
fit  low  resolution  Images  quite  well,  are  best  realized  by  elliptic  PDE 
models. 

Our  results  on  image  restoration  have  demonstrated  that  each  class  of 
PDEs  leads  to  a different  type  of  computational  algorithm.  For  example, 
the  hyperbolic  equations  yield  causal  vector-recursive,  Kalman  filtering 
type  algorithms;  the  parabolic  equations  yield  semi-causal  algorithms , i. e . , 
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algorithms  which  are  causal  recursive  in  one  dimension  and  noncausal  in  the 
other;  the  elliptic  equations  yield  noncausal  algorithms.  These  different 
structures  were  shown  to  give  new  and  different  computing  architectures  as 
well  as  performance  bounds  on  the  processing  technique. 


3. 1 Causal  Models  and  Hyperbolic  PDEs 

Consider  the  operation  equation  defined  on  (0,®)  x (0,®) 


V = 


3xSy 


+ O'.  -1^  + O'-  ~ + a.v 
1 ax  2 ay  3 


(9) 


where  v = v(x,y).  Here  is  a second  order  hyperbolic  PDE  operator  such 

that  ^|^v(x,y)  = f(x,y)  is  a well  posed  initial  value  problem.  A finite 
difference  approximation  of  this  equation  gives  a discrete  random  field 
representation  of  the  type  (to  be  called  Cl) 


“i.j 


= a.u.  , . + a-u.  . , -a.u  , < i "*■  < 

1 i-l,j  2 i,j-l  3 i-l,j-l  i,j 


(10) 


This  is  a well  known  model  used  for  DPCM  coding  of  images.  For  a^  = ^^^2 
this  equation  is  an  exact  realization  of  the  separable  covariance  model 
given  by 

R(k,£)  = . 

2 2 

where  E[u.  J = 0 and  E[u  ] = a . 

2 

The  various  parameters  of  (10),  are  identified  as  a^^  = a^  = p , = p , 

= 0,  = (1  - P^)^6k.O^A.O  ^,m  the 

Kronecker  delta  function. 
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3. 2 Semlcausal  Models  and  Parabolic  PDEs 


Consider  the  heat  equation  operator 


V + av 

XX 


(U) 


This  is  a second  order  parabolic  PDE  operator  and  leads  to  a well  posed 
initial  value  problem  in  y dimension  and  well  posed  boundary  value  problem 
in  X dimension.  A finite  discrete  approximation  of  this  gives  a semlcausal 
representation  (SCI),  which  is  causal  in  j variable  and  noncausal  in  i vari- 
able, of  the  form 


u.  = ffCu.  , , + u. .)  + Yu.  . , + 
i.,k  1-1, j 1+1, j i,J-l 


ij 


where  a < ^ |y|  < 1 and  20?  + Y < I 


(12) 


This  and  other  semlcausal  models  lead  to  computationally  desirable,  the  sd- 
called,  hybrid  algorithms.  A hybrid  algorithm  is  a recursive  algorithm  ob- 
tained after  transforming  each  row  (or  column)  of  the  image  by  a unitary 
transform.  For  such  models,  it  has  been  shown  that  a 2-dimensional  image 
can  be  decomposed  into  a sequence  of  one  dimensional,  independent  Markov 
processes.  The  statistics  of  certain  semlcausal  models  have  been  shown  to 
fit  high  resolution  images  quite  well  L‘^,6]. 


3. 3 Noncausal  Models  and  Elliptic  PDEs 

Consider  the  well  known  Poisson  operator  equation 


V + V + av 
XX  yy 
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(13) 


This  is  a boundary  value  elliptic  PDE.  Such  equations  have  been  shown  to 
yield  noncausal  stochastic  image  models  and  have  shown  to  fit,  quite  well, 
satellite  and  other  remotely  sensed  images,  which  have  nearly  isotropic 
covariance  functions.  Noncausal  models  have  also  been  useful  for  represent- 
ing images  with  separable  covariance  functions  f4,b3.  The  adavantage  of 
such  representations  is  that  often,  they  yield  computationally  fast  Kar- 
hunen  Loeve  (KL)  expansions  of  image  data  and  have  proven  useful  in  resto- 
ration and  data  compression  of  images. 


3.4  Application  to  Rdar  Image  Area  Correlation  Techniques 

For  the  sake  of  clarity  in  understanding,  we  shall  work  here  with  one 
dimensional  Images  defined  over  an  interval  Correspondingly  the 

area  correlation  is  defined  as 


c(x) 


a 

1^(1  + x)l2(5)d5 
-a*^ 


(14a) 


In  practice  the  reference  image  1^^  could  be  considered  as  representing  the 
terrain  features  very  accurately  and  a blurred  image  of  the  terrain 
given  as 

= \ n(x  - ^)I^(?)dC  , (I4b) 

'■'-a 

where  h(x)  is  the  point  spread  function  of  the  imaging  optics.  As  a first 
approximation  h(x)  is  assumed  to  be  spatially  Invariant.  Extension  to  spa- 
tially variant  case  is  possible  (although  the  computational  complexity  may 
increase).  Combining  (I4a)  and  (14b)  we  get 
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:(x)  = ^ h(^  - + §)d5d§' 


(14c) 


For  evaluation  of  an  area  correlation  system,  generally,  the  first  two  mo- 
fjients  of  c(x)  are  needed.  From  above  these  are  obtained,  assuming  Ij^  to  be 
a stationary  process  whose  autocorrelation  is  as 


i^(x)  ^ E[c(x)l  = ^a  h(5  - - l'+  x)d§d5' 


(15a) 


,(x)  = E[c^(x)>  ia  ^ h(5  - r)h(ll  - TlOE[I^(r)I^(x  + 5)I^('n  Ol^(x+Tl)l . 


dldS^dTldTl ' (15b) 


Equations  (15a)  and  (15b)  can  be  evaluated,  at  least  in  principle,  if  the 
second  and  fourth  order  statistics  of  the  random  process  are  known. 
Measurement  of  such  statistics  from  available  data  is  impractical  even  when 
stationarlty  and  ergodicity  assumptions  are  made.  This  is  because  in  (15b), 

3 

for  an  N point  discretization  of  Ij^(x),  at  least  N coefficients  for  the 
fourth  order  expectation  will  be  required.  For  two  dimensional  images, 
the  dimensionality  will  be  O(N^).  Hence,  it  is  desirable  to  have  a dif- 
ferent stochastic  representation  for  so  that  all  the  statistical  para- 
meters can  be  calculated  from  it.  One  way  is  to  determine  a stochastic  PDE 
representation  of  the  two  dimensional  image  Ij^(x,y)  and  determine  all  the 
statistics  from  that  model.  This  approach  would  be  the  2 dimensional  equi- 
valent of  calculation  df,  say,  the  autocorrelation  function  of  a one  dimen- 
sional random  process  from  its  stochastic,  ordinary  differential  equation 
model.  The  difference  would  be  that  one  would  have  to  solve  PDEs  for  two 
dimensional  variables. 


42 


IV.  RKAL-TIME  RF  SIMULATION  OF  RE-ENTRY  VEHICLES 


The  availability  of  an  image  model  is  particularly  useful  in  real  time 
simulation  of  radar  imaging  systems.  For  example,  in  a hardware  in  the  loop 
simulation  of  a terminal  guidance  system  at  the  Advanced  Simulation  Center 
of  U.S.  MIRADCOM,  Redstone  Arsenal,  Alabama,  it  is  planned  to  generate,  in 
real  time,  image  signals  as  shown  in  Fig.  28.  Here,  the  re-entry  vehicle 
(RV)  flight  hardware  (Radar  scanning  antenna,  the  radar  receiver  electronics, 
and  the  correlator  unit)  is  mounted  at  one  end  of  a chamber  and  the  scanning 
antenna  (i.e.,  the  FLR)  is  allowed  to  rotate,  as  it  operates  in  the  real 
environment.  This  radar  antenna  scans  a large  array  of  antennas  located 
at  the  other  end  of  the  RFSS  chamber.  Digital  image  signals  generated  by  a 
digital  computer  are  converted  to  analog  video  signals  and  are  put  on  Radio 
Frequency  (RF)  carrier  pulses  for  transmission  from  the  array.  By  a sophis- 
ticated computer  controlled  switching  method,  signals  are  radiated  from 
the  array  such  that  they  are  synchronous  with  the  scanning  antenna  position. 
The  signal  power  levels  and  array  position  control  are  fine  enough  to  simu- 
late the  real  time  radar  return  signals  received  by  the  RV.  The  reference 
and  the  live  FLR  Images  are  updated  in  real  time  as  the  trajectory  of  the 
RV  changes.  A sequence  of  reference  images  is  known  in  advance  at  few 
fixed  points  on  the  trajectory.  The  reference  image  at  any  point  is  prepared 
off  line  and  the  relevant  Images  stored  in  advance  for  signal  generation. 

Fig.  29  shows  a more  detailed  design  architecture  developed  by  the  prin- 
cipal investigator  during  an  LRCP  research  effort  for  a real  time  hard- 
ware in  the  loop  simulation  (at  U.S.  MIRADCCM)  of  an  area-correlation  system. 
The  principal  advantage  of  these  simulations  is  that  by  including  complex 
hardware  in  the  simulation  loop,  realistic  flight  test  data  can  be  obtained 
for  system  evaluation.  This  will  of  course  minimize  the  number  of  actual 
test  flights  needed  and  will  therefore  reduce  test  costs. 
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FIGURE  28  OVERVIEW  OF  REAL  TIME  HARDWARE  IN  THE  LOOP  SIMULATION  OF  AN  AREA  CORRELATION  GUDIANCE  SYSTEM 


FIGURE 
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ABSTRACT 

Stochastic  representation  of  discrete  Images  by 
partial  differential  equations  proposed  earlier  In  [1] 

Is  used  for  filtering  of  noisy  Images.  Comparisons 
among  various  representations  and  their  relative 
effectiveness  in  modeling  the  actual  statistics  of  the 
images  Is  shown  by  examples.  Certain  image  decompos- 
itions resulting  in  fast  filtering  algorithms  are 
given.  The  results  are  also  compared  with  spatial 
averaging  and  Fourier  domain  Wiener  filters.  Results 
show  superiority  of  noncausal  and  serolcausal  models 
over  causal  and/or  separable  covariance  models. 

I.  INTRODUCTION 

In  an  earlier  paper  (IJ  image  representation 
bv  partial  differencial  equations  (PDEs)  was  considered. 
It  was  shown  that  corresponding  to  the  PDE  classifica- 
Clon  of  hyperbolic,  parabolic , and  elliptic  systems, 
three  different  stochastic  rcpresciiLations  viz.,  causal, 
scmlcausal ,and  noncausal  representations  are  realized. 
These  representations  have  different  spatial  structures 
and  are  realizations  of  different  type  of  covariance 
functions.  For  example,  certain  elliptic  equations  can 
realize  certain  isotropic  covariance  functions  which 
cannot^®  realized  by  hyperbolic  or  parabolic  equations. 
On  the  other  hand  certain  covariance  functions  could  be 
realized  by  all  the  three  classes  of  representations. 

Use  of  PDEs  directly  for  image  representation  bypas- 
ses the  difficult  problem  of  two  dimensional  spectral 
factorization.  The  experimental  results  of  [1,2]  have 
shown  that  the  covariance  functions  realized  by  some 
of  the  PDEs  approximate  the  actual  covariance  of  dif- 
ferent image  data  better  than  the  commonly  used  covari- 
ance function  models  for  those  image  data. 

In  this  paper  we  consider  the  problem  of  linear 
filtering  of  noisy  images  represented  by  different  PDE 
models.  It  is  shown  that  the  causal,  semicausal,  and 
noncausal  models  give  rise  to  recursive,  hybrid  or 
semi-recursive,  and  transform  domain  or  nonrecurslve 
filtering  algorithms.  Recursive  filtering  for  images 
has  been  considered  earlier  by  many  researchers  13-7]. 
Typically  in  these  approaches,  one  insists  on  a causal 
image  representation  e.g.,  state  variable  or  two  dimen- 
sional autoregressive  models,  such  that  recursive  fil- 
tering equations  of  Kalman  and  Bucy,  or  equivalently, 
those  of  Aasnaes  and  Kailath  (8]  arc  applicable.  By 
hybrid  or  semirecursive  filtering  we  mean  a filtering 
method  which  is  non-recurslve  in  one  of  the  image 
demensions  and  is  recursive  in  the  other.  Algorithms 
of  this  type  include  those  vector  recursive  algorithms 
where  the  state  variable  vector  contains  an  entire 
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column  (or  row)  of  the  image,  and  have  been  considered 
by  Jain  (9],  Murphy  and  Silverman  (10]  and  others  (11, 
12).  Nonrecurslve  or  transform  domain  algorithms  arise 
quite  naturally  for  noncausal  models  and  lead  to  smooth- 
ing equations  very  similar  to  discrete  Fourier  Trans- 
form (DFT)  domain  Weiner  filters.  Similar  algorithms 
have  been  studied  previously  In  119|. 

Although  some  of  the  image  models  (and  their  assoc- 
iated filter  equations)  considered  here  have  been 
studied  earller(9,13,l9l,thls  paper  reconsideres  these 
as  well  as  some  new  models  in  the  more  general  framework 
of  PDEs.  This  has  several  advantages.  First,  the  well 
developed  theory  of  PDEs  offers  many  analytical  answers 
to  problems  related  to  stability,  stationarity , perfor- 
mance bounds,  etc.  Second,  the  different  existing  nu- 
merical methods  for  solving  PDEs  hold  promising  possi- 
bilities for  implementing  filter  solutions.  Finally, 
the  PDE  theory  lends  Insight  into  certain  filtering 
phenomena.  An  example  of  this  is  shown  in  section  b 
to  explain  the  performance  of  a spatial  averaging 
filter.  Several  filtering  experiments  on  actual  Images 
and  comparisons  of  different  models,  old  and  new, 
exemplify  the  utility  of  our  approach. 

In  section  2 we  briefly  review  the  results  of  HI 
that  are  used  here.  Specifically,  the  different  models 
identified  in  (1)  are  summarized  and  some  notation  is 
established.  In  section  3,  4 and  5 the  filtering  algor- 
ithms associated  with  the  causal,  semi-causal,  and 
noncausal  models  respectively,  are  described.  Filteilng 
experiments  on  four  different  images  have  been  performed. 
In  section  6,  the  results  of  these  experiments  and  their 
mutual  comparisions  as  well  as  with  standard  Fourier 
domain  Wiener  filtering  techniques  are  discussed.  Con- 
clusions and  possible  extensions  of  this  work  are  given 
in  section  7. 


2.  PDE  IMAGE  MOUELJ 

Stochastic  PDEs  have  been  studied  earlier  for 
various  applications  [14-18].  In  a recent  paper  Jain 
[1]  has  studied  the  possibilities  of  modeling  actual 
image  data  by  PDEs  and  many  different  models  have  been 
suggested.  Table  I sutnnarlzes  some  of  those  models 
which  we  will  be  studying  here  for  filtering  applica- 
tions.Wc  start  with  a PDE  operator  l(3/3x,3/3y)  and 
approximate  It  by  discrete  operator  D(®x  , *2K®^t.ained 
by  a suitably  chosen  finite  difference  method.  Then 
we  assume  an  image  representation  of  the  form 


where  u . represents  the  (relative)  brightness  of  the 
zero  meaA-^unlt  variance  image  at  the  point  (i,j)  for 
and  j is  a two  dimensional  zero  mean  dis- 
crete white  noise  or  some  sort  of  moving  average  process 
with  absolutely  continuous  spectral  density  function. 
Clearly,  and  in  D(e^,S2)  also  be  chosen  to  be 
the  two  dimensional  ^ - transform  variables.  If  we  de- 
note spectral  density  function  of 
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(c.  .}  for  •.-exp(Jw  ),  • - exp  (Jw^),  then  the  spec- 

tral’^denslty^function  of  ' 


j ) Is 


tlven  by 


S (o, 


.) 


1 


i covariance  function  of  {u^  ) is  obtalneH  quite 

sily  by  the  Fourier  inverse  ol^S  (e-^'^1,  e^'*2),  which 
could  be  accomi  wished  by  a suitable  FFT  algorithm. 

Table  I lists  the  forms  of  chosen  for  various  models. 
The  values  of  different  parameters  identified  for  these 
models,  for  a 236x236  Girl  image, are  listed  in  Table 
II.  Here  represents  K(c  ] and  is  the  value  of 
t f\  f\\  T ..U...  .1  » J y. 


V..  ..w.w  ..  ^ i 

(0,0).  Table  I also  shows  *tlhe  spatial  structure  of 
the  different  models  arising  from  a finite  difference 
approximation  of  the  operator  1.  Here,  if  u corres~ 
ponds  to  the  point  marked  A,  then  the  spatial 'structure 
shows  all  other  points  that  are  used  directly  in  the 
prediction  of  A.  The  stability  of  various  models  and 
their  stationarity  (asymptotically)  is  assured  by  im~ 
posing  restrictions  on  the  model  parameters  which  are 
consistent  with  the  well  posedness  of  the  original 
PDEs.  The  parameters  of  Table  II  satisfy  those  re-* 
strictions  and  realize  stationary  fields.  Parameters 
for  some  other  images  and  related  details  can  be  found 
in  (1,21. 

We  will  consider  the  simplest  filtering  problem 
where  the  observations  are  given  as 


'i.j 


i.j 


'’i.J* 


(1) 


CAUSAL  MODELS 


Causal  aodela  considered  here  result  from  a 
hyperbolic  system  of  PDEs  and  have  spatial  structure 
of  the  form  ' 


i.J 


Vi-i,j 


^'•2“i.j-rVi-i.j-i^i.j* 


(2) 


For  a^  - aj^a2  this  results  in  an  often  used  separable 


covariance  function.  Such  a model  has  been  con- 
sidered by  Hablbi  |6)  for  recursive  filtering)  but  the 
equations  developed  there  have  been  shown  to  be  in 
error  by  Strintzis  (7J.  Recently  Bary,  et  al  [17)  have 
developed  some  asymptotic  recursive  filters  for  such 
models.  More  recently  Wong  116]  has  mentioned  that  two 
dimensional  recursive  filter  for  the  continuous  equiva- 
lent of  this  model  exists  and  that  the  filtering  equa- 
tions are  rather  complicated.  Here  we  consider  the 
recursive  filtering  equations  for  the  discrete  model 
of  (2)  and  show  that  theaa  equations  may  be  obtained 


from  a set  of  one  dimensional  vector  recursive  equat- 
ions . 

If  u denotes  the  Nxl  column  vector  of  an  NxM  image 
then  (2)  can  be  written  in  vector  form  as 


(3) 


where  and  are  nonsingular  lower  triangular  Toe- 
plltz  matrices  and  b are  Nxl  vectors  containing  initial 
conditions.  If  b is  orthogonal  to  u . and  t,  and  if 
it  is  a white  nolde  process  (for  a -ajS*  this  is  true) 
and  if  we  define  f.  ■ t.  + b then'*(3)  becomes  a vector 

J J 4 


Markov  process 


i 


(4) 


The  observations  of  (I)  can  be  written  in  a vector  form 
as 


Yj  ■ u^  + Oj 


(5) 


Various  vector  recursive  filtering  equations  can  now  be 
written  quite  straightforwardly.  For  example,  if  a 
column  vector  is  scanned  at  a time  and  the  previously 
scanned  vectors  are  used  to  predict  the  next  column 
vector,  then  the  Kalman  filtering  equation  for  one 
step  predictor  can  be  written  as 


.-1 


■^j+l*  ‘■I  ‘■2  V° 


(6) 


where  (i,j)  belongs  to  a set  I,  of  integer  pairs 
and  (n.  .}  2^^  ^ zero  mean  white  Gaussian  noise  of 
variance^  . The  set  I will  be  suitably  defined  for 
different  Image  models  as  we  proceed  in  the  paper.  We 
would  be  Interested  in  the  various  predictors,  on-line 
filters  and  smoothers  of  u , given  the  observations 
of  (1).  It  should  be  noted'that  the  models  are 
stationary  for  Infinite  images.  To  apply  them  for  re- 
storation of  finite  size  images  certain  Initial  and/or 
boundary  conditions  on  the  models  will  have  to  be 
specified.  If  stationarity  is  to  be  retained,  the 
boundary  conditions  should  be  consistent  with  the  cov- 
ariance function  realized  by  the  model.  In  practice, 
these  consistent  specifications  may  be  difficult  at 
times  and  one  may,  instead,  settle  for  a nonstationary 
model  which  is  asymptotically  stationary,  by  assuming 
a more  convenient  model  for  the  boundary  variables. 

In  fact  we  will  see  that,  implementation  wise,  some  of 
the  simplest  algorithms  are  obtained  for  nonstationary 
(but  spatially  Invariant)  models  whose  boundary  con- 
ditions are  deterministic.  This  would  result  in  sub- 
optimality,  which  depends  on  image  size,  and  may  not 
always  be  significant  compared  to  the  error  in  modeling. 
Throughout  the  paper  we  will  attempt  to  point  out  such 
optimality  vs  feasibility  tradeoffs. 


where  G.  is  Kalman  gain  matrix  that  is  obtained  by  ■ 

solving*^ an  associated  Riccati  equation.  For  a true  < 


two  dimensional  recursive  filter  suppose  the  image  is  a j 
scalar  scanned  image,  e.g.,  raster  scanning,  say  from 
top  to  bottom  (along  i)  and  left  to  right  (along  j). 

For  such  a scanning  the  one  step  predictor,  u 
is  defined  as  the  best  linear  estimate  of 
the  observations  upto  the  previous  scan  point  '^''^i.e., 
y . The  solution  of  such  a recursive  filter,  results 
ln*a  one  step  predictor  given  by 


i+l,j 


Xj(l+0 


<7) 


where  Xj  Is  a "state  vector"  for  the  estlnates  at 


8canner'’po8ltlon  (l,j)  and  satisfies 


1+1 


x“  - Uj,  OslxN-1  (8) 


1-1 


I.J  ^i.J 


- Xj  (1) 


(9) 


Here  x.  la  a Nxl  column  vector  representing  the  best 
estimate  of  jth  column  vector  u given  observations  only 
up  to  ; v^  la  the  scalar,  scanner  Innovations  1 

process, and  fal  Is  another  gain  vector  and  Is  a fun- 
ction of  G of  (6)  Note  that  the  Initial  condition  of 
(8)  requires  first  solving  the  vector  recursive  filter 
equation  (6)  up  to  step  J.  Hence,  we  have  to  solve  two 
sets  of  vector  equations  l.e,  , (6)  and  (8)  one  In  the 
horizontal  (j)  direction  and  the  other  In  the  vertical 
(1)  direction  to  obtain  a two  dimensional  scanner  re- 
cursive filter.  The  major  computational  effort  Is  In 
the  solution  of  the  Rlccatl  equation  for  C and  h^.  For 
large  Image  size  G.  and  h^  can  be  replaced  by  their 
steady  state  values'^and  the  resulting  filter  equations 
then  require  about  O(N^)  operations  per  vector  step 
(or  a total  of  about  O(N^)  operations).  Further  simpli- 
fications are  possible  by  noting  that  G becomes  a 
Toeplltz  matrix  when  1,J  and  (6)  rednces  to  convol- 
ation  operations.  For  the  separable  model  (l.e.  for 
a3  - ■ia2)  en  Iterative,  smoothing  filter  Is  also 
possible.  If  we  write  (1)  In  matrix  form 


U - Y + n. 


(10) 


and  if  j.  ^ than  the 


smoothed  estimate  of  the  Image  using  Cl  model  can  he 
found  bv  iterative  solution  of  the  following  equation 
In  the  sine  transform  domain  [2J* 

2 


Q u. 


If  - B,..  U 

k o 


-g-t 

2, 

a”P  / (1-p  ) , 


(DU.  Q4<JU.  D) . 


and  Q and  D are  NxN  matrices  defined  by 


li-Jl  -1 


QUj-PUj.i+tj+bj, 


li  JS  N 


where  Q is  an  HxH  synunetrlc,  tridiagonal,  Toeplltz, 
matrix  defined  In  (13);  P « yl  for  SCI  and  P^pQ  for 
SC2  models;  and  b.  is  a Nxl  colmn  vector  consisting  of 
only  boundary  valoes  identified  as 


'a.(Uo  j,0. 


, j)  » for  SC1« 


^“o,J-l’° “n+1,j‘‘’“n+1,J-i\  for 


The  statistics  of  could  then  be  identified  as 


Etc  e^l 


. for  SCI  & SC2, 

b'i  «j.k- 


Q , for  SC2 

J 


The  semlcausal  sturcture  of  (lA)  Is  evident  by 
noting  the  noncausallCy  of  any  element  u.  . of  the 
linage  In  the  i variable  and  its  causality’in  J varia- 
ble. Any  attempt  to  reindcx  these  equations  as  re- 
cursive in  the  1 variable  would  result  in  unstable 
causal  systems  and  will  no  longer  represent  (asympto- 
tically) stationary  random  fields  [9]. 

Letting  b contain  arbitrary  boundary  values, 
randog  or  detenlnlstlc,  we  define  two  processes  u 
and  Uj  as  ^ 


“j  - “j.i  ^ -y 
Q uj  - P uj.^  . b^. 


u°.0 


For  a nonsingular  matrix  (),  (lA)  will  have  a unique 
solution.  From  this  uniqueness  of  u.  and  linearity  of 
(14)  It  follows  from  (16)  and  (17)  that  u.  has  the 
stochastic  decomposition  ^ 


otherwise. 


d,  , ■ d,,  „ • 1,  d.  , ■ 0,  otherwise,  where 
1,1  N,N  ifj 

2 

o “ p/(l+p  ).  The  reason  we  opt  to  solve  (11)  and  (12) 
In  sine  transform  domain  is  that  the  sine  transform 
diagonallses  the  matrices  Q and  D and  then  they  could 
be  solved  as  set  of  N scaler  Independent  equations, 
thus  resulting  In  computational  savings.  Details  of 
derivation^  will  be  given  elsewhere. 


4.  SEMICAUSAL  MODELS 

If  u denotes  the  Nxl  column  vector  of  an  NxN 
image  then,  the  semlcausal  models  SCI  and  SC2  given  in 
Table  1 can  be  written  as 


Uj  - u°  + Uj  , 1 * j * N.  (18) 

The  vector  u^  is  called  the  boundary  response  of  u and 
is  completely  determined  by  the  initial  and  the  boondary 
values  of  the  image.  This  decomposition  would  be  orth- 
ogonal if  (bj}and{c^}  are  uncorrelated  sequences  i.c., 

■ “ "j... 

The  above  condition  is  always  satisfied  for  the  SC2  model. 
Note  that  (16)  is  now  a nonstationary  Harkov  process 
starting  with  a zero  initial  state.  Using  (18)  in  eqn. 

(5)  the  observations  can  be  written  as 

y°  - Pj  + uj  (19) 


y"  - “j  + • (20) 

Now  noting  the  fact  that  the  discrete  sine  transform 
Y diagonallses  any  tridiagonal  symmetric  Toeplltz  mat- 
rix, such  as  P andQ/he  vector  equations  (16),  (17),  (19) 
and  (20)  can  be  converted  into  scaler  decoupled  set  of 
causal  Markov  representations  by  taking  their  sine 
transform.  The  Kalman  filtering  equations  for  the  re- 
sulting equations  in  the  Sine  transform  domain  are 
given  quite  easily  as  follows. 

One  Step  predictor 


■(1  + Yj/P,^^)  ^ 


ElCj  (1)1 


On-line  filter 


(e.  - V, 


Interpolator  (Smoother) 

V* 


c.fj.Wj^^  «j(.j-Vj)/a^ 


Yyj  . Vj  fuj,  Cj 

c(l)  - y/1^  and  d(l)  - f 


c(i)  - 0 


and  d(l)  * A. 


for  SCI, 
for  SC2. 


All  the  above  equations  for  filters  are  decoupled  in 
the  variable  i and  can  be  solved  independently  for 
each  i.(The  subscript  i has  been  omitted  for  notational 
simplicity).  Prom  these  the  spatial  domain  image  es- 
timates are  given  by  Inverse  unitary  transformation 
(fl  • . f),  i.e. 


.o*.  b 

U +Uj 


* b 

TVj  + ul“  , etc. 


The  semlcausal  model  recursive  algorithm  Is,  very  sim- 
ply, as  followa.  From  the  given  boundary  values  (or 


•The  sine  transform  Is  related  to  the  discrete  Fourier  transform  and  has  a fast  Implementation  algorithm  similar 
to  FFT  algorithm.  siii.a  » 
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Citeir  estimates)  first  calculate  u \ and  determine  the 
nkHllflod  observations  y . Then  transform  y by  the 
(lirsl)  Sine  transform  and  solve  the  desirea  recursive 
filter  equations.  The  inverse  transformation  of  the 
recursive  filter  estimate  and  adding  u^  gives  the  de- 
sired filtered  image  vector.  For  details  of  SC2  model 
see  19]. 


tic.  In  that  case*  these  are  estimated  first  and  are 
used  to  determine  This  will  make  our  scheme  slj- 

ghtly  subopt imal  although  its  signilicance  decreases 
as  the  image  size  increases.  See  |19]  for  an  analysis 
of  boundary  effects. 

6.  EXAMPLES  AMD  DISCUSSION 


5.  NONCAUSAL  MODELS 


From  the  spatial  structures  of  the  noncausal 
models  NCI,  NC2  and  NC3  (see  Table  1)  It  can  be  seen 
that  they  are  nonrecursive  in  both  the  indices  I and  j. 
Any  attempt  to  reindex  them  so  as  to  yield  a recursive 
relation  in  either  of  the  indices  would  result  in  an 
unstable  system.  Let  U be  an  NxN  Image  matrix  and  1) 
be  an  N^xl  vector  obtained  by  the  lexicographic  order- 
ing of  U,  then  the  noncausal  models  can  be  written  In 
Che  form 


J U 


(24) 


2 2 - 2 

where  J is  an  N xN  matrix*  c is  the  N xl  vector  cor- 
responding to  Che  matrix  t containing  and  B 

Is  a N^xl  vector  whose  entries  depend  only’in  the 
boundary  values.  From  (24)  one  gets  Che  decomposition 


U 


U°  + u'’  or  U 


where 


u + u 


j"^  B. 


(25) 

(26) 


r.b 


Note  that  U is  completely  determined  from  the  boundary 
values  and  Is  called  the  boundary  response.  Let 


A F [e.e’’^]  - 


^ . 


then  from  (26)  and  (27)  the  covariance  matrix  of  U 
obtained  as  . 


R°A  ElU°.U'’ 


(27) 
is 

(28) 

The  structures  of  the  matrices  J and  t resulting  from 
the  PDFs  considered  is  such  that  the  Kronecker  product 
the  Sine  transform  y contains  all  the  eigen- 
vectors of  J and  I for  NCI  and  NC3  models  (For  NC2 
model  this  is  still  a good  approximation).  Therefore 
yffl'l'  becomes  tlie  Karhunen  Loeve  (KL)  transform  of  D° . 
Hence*  if 

V°  - U°  or  V°-  (29) 

then  V°  is  the  K-L  transformations  of  U°  and  its  ele- 
ments will  be  uncorrelated.  As  before,  the  observa- 
tions can  also  be  decomposed  to  give 


where 


Y + U 


u + n 


(30) 

(31) 


Now  if  the  boundary  values  (or  equivalentlv  U^)  were 
deterministic  we  can  obtain  the  Wiener  filter  equa- 
tions for  (23)  and  (31)  in  KL  transform  (Y)  domain  as 
N^  decoupled  scaler  equations  and  are  found  to  be 


l.j 


6^+0  /k 

n 1,J  i.j 


Isi.jsN. 

2. 


where  V • fUf,  2 - m,  C ■ fBf  , o - E[n.  ].< 
and  k are  the  eigenvalues  of  J and  ^ 

I respectively.  The  optimal  smooched  estimate  of  j^he 
Image  is  then  found  by  Inverse  sine  transform  of  V , 
l.e.  * 


U 


tv  T. 


In  practice  Che  boundary  values  may  not  be  determinls- 


The  filtering  techniques  discussed  In  the  pre- 
vious sections  were  implemented  on  four  65x63  blocks 
of  four  236x236  Images  corrupted  by  white  Gaussian 
noise.  The  Signal  to  Noise  ratio  (S/N)  values  of  1* 
2*  and  3 were  used.  This  ratio  is  defined  as 

S/fj  a Standard  deviation  of  the  signal 
Standard  deviation  of  the  noise 


The  various  models  are  compared  mutually  as  well  as 
with  their  corresponding  Fourier  domain  Wiener  filters 
and  spatial  averaging  filters.  The  spatial  averaging 
filter  averages  each  pixel  with  the  nearest  four  points* 


l.e., 


i.j 


For  causal  model  Cl  the  Iterative  solution  of  (11) 
and  (12)  was  used.  For  solving  semicausal  and  non- 
causal filter  equations  the  boundary  values  were  first 
estimated  separately.  For  these  models  the  size  of 
the  image  block  U was  taken  to  be  63x63  and  with  boun- 
daries the  total  block  size  becomes  63x63.  The  four 
boundary  vectors  were  estimated  independently  using  one 
dimensional  Fourier-Wiener  filter  and  first  order  Markov 
scati.stics . 

These  boundary  estimates  were  used  (instead  of  the 
actual  boundary  values)  In  computing  the  boundary  re- 
sponse functions..  For  NC2  model  two  boundaries  on  each 
side  of  Che  image  are  needed  and  we  have  assumed  the 
outer  most  boundaries  to  be  zero. 

Tables  HI*  IV*  and  V show  the  comparisons  among 
Che  performance  of  various  filters  for  S/N  of  1,2  and 
5 respectively.  The  Cl,  SC2  and  NC3,  all  correspond 
to  the  separable  covariance  model  and  hence  their  re- 
sults differ  only  In  terms  of  the  requirements  of 
knowledge  of  the  boundaries.  As  expected,  their  per- 
formance Is  quite  close  and  the  difference  depends  on 
how  good  are  the  boundary  estimates.  Their  poor  per- 
formance for  the  blocks  of  the  MOON  and  ERTS  images 
was  found  to  be  due  to  the  fact  that  the  one  step  cor- 
relation parameter,  p,  for  these  blocks  was  significan- 
tly different  from  the  one  measured  for  the  larger 
256x236  Images.  On  the  other  hand,  the  other  models 
(SCI,  NCI  and  NC2)  perform  relatively  very  well  on 
these  blocks,  This  implies  that  the  latter  models 
are  much  more  stable,  l.e.,  they  are  relatively  In- 
sensitive to  the  local  changes  in  statistics. 

The  superior  performance  of  the  NCI  and  SCI  models 
is  evident  from  their  filter  performances.  Except  for 
the  couple  image  block  the  performance  of  NCI  is  the 
best  in  almost  all  the  cases.  In  most  cases  NC2  per- 
forms better  than  the  separable  models  but  is  poorer 
than  NCI  and  SCI  models.  The  superior  performance  of 
these  models  may  be  attributed  to  the  superiority  of 
their  fitting  to  the  image  statistics. 

Tlie  simple,  adhoc.  Spatial  Averaging  filter  performs 
surprisingly  well,  especially  for  moderate  signal  to 
noise  ratios.  Note  that  it  performs  better  than  the 
Cl,  SC2  and  NC3  models.  This  can  be  explained  by  th« 
fact  that  the  transfer  function  of  the  spatial  averag- 
ing filter  Is  a fairly  good  approximation  of  the 
Fourier-Wiener  filter  corresponding  to  the  NCI  model, 
which  a.s  we  have  seen  here.  Is  a very  good  model  for 
filtering  of  images.  Tables  III  to  V also  show  that 
almost  ail  the  filters  perform  somewhat  better  than 
their  corresponding  Fourier  domain  Wiener  filters. 

This  Is  because  for  finite  size  Images,  the  Fourier 
domain  samples  (of  observations  and  the  original  data) 
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are  correlateJ  and  the  corresponding  Wiener  filter  Is 
subopclmal  120].  The  aeolcausal  and  noncousal  filter 
are  also  suboptlmal  whenever  the  boundary  values  are 
not  known  exactly.  However,  this  suboptlmallty  Is  less 
significant  than  the  Fourier  domain  suboptlmallty. 

Fig.  1 shows  some  of  the  Images  resulting  from 
the  application  of  the  foregoing  filters.  The  perfor- 
mance differences  due  to  different  filters  were  also 
evident  on  a visual  display. 


7.  CONCLUSIONS 


In  conclusion.  It  was  demonstrated  that  different 
FDE  classes  yield  different  filter  structures  and  also 
different  filter  performances.  Many  of  these  models 
offer  substantial  advantages,  both  In  terms  of  compu- 
tational complexity  as  well  as  filter  performance,  over 
the  conventional,  separable  Image  covariance  models. 

More  general  PDE  models  may  be  used  to  achieve  even 
better  results.  Mso,  these  models  may  be  used  for 
other  two  dimensional  signal  processing  applications 
such  as  data  compression,  detection,  synthesis, etc.  [21] . 
Inevitably,  In  all  these  applications  numerical  methods 
for  solving  PDEs  will  play  a major  role. 
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a)  Original 


b)  Noisy 


c)  Separable  Model  F-W  Filter 


d)  NCI  Model  F-W  Filter 


e)  SCI  Model  Filter 


f)  NCI  Model  Filter 


Fig.  1:  SMOOTHING  OF  A 64  x 64  BLOCK  OF  THE  GIRL  IMAGE 
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SUMMARY  OF  SOME  PDE  MODELS  AND  THEIR  DISCRETE  APPROXIMATIONS  D[U]  - c USED  FOR  IMAGE  MODELLINC. 


Covariance  same  as 





f 


TABLE  II 

ACTUAL  PAKAMETERS  MEASURED  FOR  THE  2S6x236  GIRL  IMAGE  TO  MATCH  THE 
16x16  COVARIANCE  MATRIX  OP  THE  IMAGE  TO  THAT  OF  THE  MODEL  FOR  VARIOUS 
MODELS. 


s. 

NO. 

MODEL 

PARAMETERS 

■ 

Causal 

CL 

p - 0.962,  . 8^  - (1-p^)^-  0.0055. 

2. 

Scmlcausal 

SCI 

a - 0.4275,  T • 0.1415,  8^-0.0198. 

■ 

Semlcausai 

SC2 

p - 0.962,  B^-(l-p^)^/l-M)^)“0.0029, 

o • p/{l-tp2) . 

■ 

Noncausal 

NCI 

a - 0.2496,0^^-0.95.  a , 8^-0.0744. 

5. 

Noncausal 

NC2 

o - 0.2494  , 8^  - 0.0053. 

■ 

Noncausal 

NC3 

p - 0.962,  B^-(l-p^)^/(l+p^)^-0.0015, 
a - p/(l+p^). 



TABLE  III 

COMPARISION  OF  VARIOUS  RESTORATION  FILTERS  IMPLEMENTED  ON  65x65  BLOCKS  OF  j 
DIFFERENT  NOISY  IMAGES  FOR  SIGNAL  TO  NOISE  RATIO  -1.  THE  ENTRIES  SHOW  THE  j 
IMPROVEMENT  IN  SIGNAL  TO  NOISE  RATIO  IN  DECIBELS. 


FILTERING 

MODEL 

GIRL 

IMAGE 

MOON 

IMAGE 

COUPLE 

IMAGE 

ERTS 

IMAGE 

Four  Point 

Spatial  Averaging 

5.23  db 

5.13  db 

5.14  db 

5.31  db 

Causal  Cl 

Four  ler**  Weiner 

9.12  db 

7.20  db 

8.10  db 

8.42  db 

(F-W)  Filter 

8.38  db 

7.32  db 

7.98  db 

8.24  db 

Sealcausal  SCI 

9.46  db 

8.89  db 

9.30  db 

9.64  db 

F-H  Filter 

8.35  db 

8.78  db 

8.32  db 

8.82  db 

Sealcausal  SC2 

8.87  db 

6.61  db 

8.27  db 

8.54  db 

F-W  Filter 

8.38  db 

7.32  db 

7.98  db 

8.24  db 

Noncausal  NCI 

1U.05  db 

8.75  db 

8.70  db 

9.88  db 

F-W  Filter 

9.00  db 

8.58  db 

7.78  db 

8.80  db 

Noncausal  NC2 

9.48  db 

8.05  db 

8.04  db 

9.32  db 

F-W  Filter 

8.24  db 

7.87  db 

7.02  db 

8.08  db 

Noncausal  NC3 

8.92  db 

6.59  db 

8.31  db 

8.46  db 

F-W  Filter 

8.38  db 

7.32  db 

7.98  db 

8.24  db 

COMPARISION  OF  VARIOUS  IMAGE  RESTORATION  FILTERS  IMPLEMENTED  OH  65x65  BLOCKS 
OF  DIFFERENT  NOISY  IMAGES  FOR  SIGNAL  TO  NOISE  RATIO  -2.  THE  ENTRIES  SHOW  THE 
IMPROVEMENT  IN  SIGNAL  TO  NOISE  RATIO  IN  DECIBELS. 


FILTERING 

MODEL 

GIRL 

IMAGE 

Four  Point 

Spatial  Averaging 

4.75  db 

Causal  Cl 

5.23  db 

Fourier- Weiner 
(F-W)  Filter 

4.59  db 

Semicausal  SCI 

5.68  db 

F-W  Filter 

4.82  db 

Seroicausill  SC2 

5.01  db 

F-W  Filter 

4.59  db 

Concausal  NCI 

6.64  db 

F-W  Filter 

5,79  db 

Noncausal  NC2 

5.77  db 

F-W  Filter 

4.42  db 

Noncausal  NC3 

5.02  db 

F-W  Filter 

4.59  db 

MOON 

1 COUIM.E 

ERTS 

IMAGE 

1 IMAGE  1 

IMAGE 

3.00  db 

3.01  db 


5.00  db 
6.82  db 


2.51  db 
3.01  db 


5.29  db 
5.15  db 


4.00  db 
3.70  db 


2.49  db 
3.01  db 


4.24  db 
3.21  db 


4.44  db 
4.34  db 


3.82  db 
3.62  db 


COMPARISION  OF  VARIOUS  IMAGE  RESTORATION  FILTERS  IMPLEMENTED  ON  65x65  BLOCKS 
OF  DIFFERENT  NOISY  IMAGES  FOR  SIGNAL  TO  NOISE  RATIO  -5.  THE  ENTRIES  SHOW  THE 
IMPROVEMENT  IN  SIGNAL  TO  NOISE  RATIO  IN  DECIBELS. 


FILTERING 

MOON 

COUPLE 

MODEL 

IMAGE 

IMAGE 

Four  Point 

3.29  db 

1.72  db 

2.80  db 

Spatial  Averaging 

Causal  Cl 

0.91  db 

-1.78  db 

0.18  db 

Fourier-Weiner 

(F-W)  Filter 

0.49  db 

-1.86  db 

0.22  db 

Semicausal  SCI 

2.11  db 

1.08  db 

2.81  db 

F-W  Filter 

1.31  db 

0.94  db 

1.80  db 

Semicausal  SC2 

0.78  db 

-1.92  db 

0.22  db 

F-W  Filter 

0.49  db 

-1.86  db 

0.22  db 

Noncausal  NCI 

2.66  db 

2.02  db 

2.24  db 

F-W  Filter 

2.54  db 

2.04  db 

2.02  db 

Noncausal  NC2 

2,06  db 

-0.39  db 

0.67  db 

F-W  Filter 

0.36  db 

-0.64  db 

-0.68  db 

Noncausal  NC3 

0.79  db 

-2.01  db 

0.22  db 

F-W  Filter 

0.49  db 

-1.86  db 

0.22  db 

1.14  db 
0.72  db 


-1.04  db 
-1.22  db 


0.08  db 
-1.09  db 


-1.02  db 
-1.22  db 


APPENDIX  II 


1.  AN  OPERATOR  FACTORIZATION  METHOD  FOR  RESTORATION  OF  BLURRED  IMAGES 

A problem  of  restoration  of  Images  blurred  by  space  Invariant  point 
spread  functions  (SIPSF)  is  considered.  The  SIPSF  operator  is  factorized 
as  a sum  of  two  matrices.  The  first  term  is  a polynomial  of  a noncirculant 
operator  P and  the  second  term  is  a Hankel  matrix  which  affects  only  the 
boundary  observations.  The  image  covariance  matrix  is  also  factorized  into 
two  terms;  the  covariance  of  the  first  term  is  a polynomial  in  P and  the  se- 
cond term  depends  on  the  boundary  values  of  the  image.  Thus,  by  modifying 
the  image  matrix  by  its  boundary  terms  and  the  observations  by  the  boundary 
observations  it  is  shown  the  Wiener  filter  equation  is  a function  of  the 
operator  P,  and  can  be  solved  exactly  via  the  eigenvector  expansion  of  P. 

The  eigenvectors  of  the  noncirculant  matrix  P are  a set  of  orthronormal 
harmonic  sinusoids  called  the  Sine  transform  and  the  eigenvector  expansion 
of  the  Wiener  filter  equation  can  be  numerically  achieved  via  a fast  Sine 
transform  algorithm  kdilch  is  related  to  the  f.ist  Fourier  transform  algorithm. 
The  factorization  therefore  provides  a fast  Wiener  restoration  scheme  for 
images  and  other  random  processes.  Examples  on  255  x 255  images  are  given. 


KEY  WORDS:  Wiener  Filtering,  Karhunen  Loeve  Transform,  Image  Restoration, 
Image  Processing 


2.  PARTIAL  DIFFERENTIAL  EQUATIONS  AND  FINITE  DIFFERENCE  METHODS  IN  IMAGE  PROCESSING 

PART  II:  IMAGE  RESTORATION 

Application  of  Partial  Differential  Equation  (PDE)  models  for  restora- 
tion of  noisy  images  is  considered.  The  hyperbolic,  parabolic,  and  elliptic 
classes  of  PDEs  yield  recursive,  semirecurslve,  and  nonrecursive  filtering 
algorithms.  The  two  dimensional  recursive  filter  is  equivalent  to  solving 
two  sets  of  filtering  equations,  one  along  the  horizontal  direction  and 
other  along  the  vertical  direction.  The  semirecursive  filter  can  be  imple- 
mented by  first  transforming  the  image  data  along  one  of  its  dimensions,  say 
column,  and  then  recursive  filtering  along  each  row  independently.  The  non- 
recursive filter  leads  to  Fourier  domain  Wiener  filtering  type  transform  do- 
main algorithm.  Comparisons  of  the.  different  PDE  model  filters  are  made  by 
implementing  them  on  actual  image  data.  Performances  of  these  filters  are 
also  compared  with  Fourier  Wiener  filtering  and  spatial  averaging  methods. 
Performance  bounds  based  on  PDE  model  theory  are  calculated  and  implementa- 
tion trade-offs  of  different  algorithms  are  discussed. 

KEY  WORDS:  Recursive  Filtering,  Two  Dimensional  Filtering,  Kalman  Filtering, 
Wiener  Filtering,  Image  Restoration,  Partial  Differential  Equa- 
tions 
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3.  FRAME  TO  FRAME  RESTORATION  OF  DIFFUSION  IMAGES 


Frame  to  frame  Image  data  is  acquired  in  many  applications  such  as  Radar, 
Biomedical,  and  Television  imaging.  In  many  situations,  the  imaging  pheno- 
menon can  be  modelled  by  a diffusion  process.  The  sequence  of  image  frames 
obtained  may  represent  the  observations  of  that  process.  Here,  we  consider 
the  problem  of  recursive  filtering  of  such  images  on  a frame  to  frame  basis. 
The  major  difficulties  are  computational  because  of  three  dimensional  struc- 
ture of  the  problem.  In  this  paper  a constant  coefficient  diffusion  system 
is  used  for  an  interframe  image  restoration  problem.  Also,  a two  dimensional 
problem  of  restoration  of  blurred  images  can  be  solved  by  imbedding  it  in  a 
three  dimensional  recursive  filtering  problem  without  blur.  The  model  struc- 
ture leads  to  a computationally  feasible  filtering  algorithm  achieving  large 
reduction  of  dimensionality  and  is  useful  in  real  time  hardware  simulation 
of  generation  of  such  blurred  image  data  as  might  occur  in  a forward  looking 
radar  (FLR). 

KEY  WORDS:  Image  Restoration,  Distributed  Parameter  Systems,  Recursive 
Filtering,  Diffusion  Equation,  Image  Simulation 


4.  FAST  INVERTION  OF  BANDED  TOEPLITZ  MATRICES  BY  CIRCULAR  DECOMPOSITION 


Banded  Toeplitz  matrices  of  large  size  occur  in  many  practical  problems 
^1-63.  Here  the  problem  of  inversion  as  well  as  solving  simultaneous  equa- 
tions of  the  type  Hx  = y,  when  H is  a large  banded  Toeplitz  matrix,  is 
considered.  It  is  shown  via  certain  circular  decompositions  of2H  that  such 
equations  may  be  exactly  solved  in  0(Nlog2N)  rather  than  in  0(N  ) computa- 
tions as  in  Levinson-Trench  algorithms.  Further,  the  algorithms  of  this 
paper  are  non-recursive  (as  compared  to  the  Levinson-Trench  algorithms)  and 
afford  parallel  processor  architectures  and  others  such  as  transversal  fil- 
ters [173,  where  the  computation  time  becomes  proportional  to  N rather  than 
NlogN.  Finally,  a principle  of  matrix  decomposition  for  fast  inversion  of 
matrices  is  introduced  as  a generalization  of  the  philosophy  of  this  paper. 

KEY  WORDS:  Toeplitz  Matrix  Inversion,  Fast  Matrix  Inversion,  Circulant 
Matrices 
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